MODULE ModArray

!!$ ModArray contains functions and subroutines that manipulate arrays.
!!$ Program conventions are from Nyhoff and Leestma, _Fortran 90 for Scientists and Engineers_, 1997
!!$
!!$ STH: Samuel Harrold

  IMPLICIT NONE

CONTAINS



!!$ Subroutine takes in 3D array with the upper and lower bounds of indices for each dimension and writes out a 2D array by mapping the 3rd dimension into the 2nd. The output is to be plotted with IDL.
!!$ Code is most versatile among compilers by declaring array size. Using an assumed array (e.g. Array3D(:, :, :) can cause problems with specialized compilers and can lead to out-of-bound references.
!!$ Note: Code assumes start of indices to be 1.

  SUBROUTINE Write3DArrayAs2DArray(Array3DIn, NumEltsD1In, NumEltsD2In, NumEltsD3In, FileNameOut)

    IMPLICIT NONE

!!$ STH: Input variables   
    INTEGER, INTENT(IN) :: NumEltsD1In, NumEltsD2In, NumEltsD3In
    REAL, DIMENSION(NumEltsD1In, NumEltsD2In, NumEltsD3In), INTENT(IN) :: Array3DIn
    CHARACTER(LEN = *), INTENT(IN) :: FileNameOut
    
!!$ STH: Internal variables
    INTEGER :: OpenStat, CloseStat, AllocStat
    INTEGER :: Row, Col, NumRowsOut, NumColsOut
    REAL, DIMENSION(:, :), ALLOCATABLE :: Array2DOut
    CHARACTER(LEN = 32) :: NumColsStr, FmtOut

    OPEN (UNIT = 10, FILE = FileNameOut, STATUS = "REPLACE", ACCESS = "SEQUENTIAL", &
         ACTION = "WRITE", BLANK = "NULL", DELIM = "NONE", FORM = "FORMATTED", &
         PAD = "YES", POSITION = "APPEND", IOSTAT = OpenStat)
    IF (OpenStat > 0) STOP "*** Cannot open file ***"

!!$ STH: Number of elements in each dimension for the output array.
    NumRowsOut = NumEltsD1In
    NumColsOut = NumEltsD2In*NumEltsD3In

    ALLOCATE(Array2DOut(NumRowsOut, NumColsOut), STAT = AllocStat)
    IF (AllocStat /= 0) STOP "*** Not enough memory ***"

!!$ STH: Array2DOut(i, j + nj*(k - 1)) = Array3DIn(i, j, k), 
!!$ where nj = number of elements in dimension 2
!!$ The starting indices for elements in Array2DOut will be 1 
!!$ regardless of the indexing scheme of Array3DIn
    Array2DOut = RESHAPE(SOURCE = Array3DIn, SHAPE = (/NumRowsOut, NumColsOut/), &
         PAD = (/0.0/), ORDER = (/1, 2/))

!!$ STH: Format output in scientific notation; 
!!$ ESw.dEe, w >= d+e+5, w = 10, d = 3, e = 2; e.g. -1.234E+56
!!$ Convert NumColsOut to a string. The maximum number of columns
!!$ that may be written out is 10^16 - 1 with FMT = '(I16.1)'
    WRITE (UNIT = NumColsStr, FMT = '(I16.1)') NumColsOut
!!$ Repeat write format for all columns in the output row.
    FmtOut = '(' // TRIM(ADJUSTL(NumColsStr)) // '(3X, ES10.3E2))'
    DO Row = 1, NumRowsOut, 1
       WRITE (UNIT = 10, FMT = FmtOut, ADVANCE = "YES") (Array2DOut(Row, Col), &
            Col = 1, NumColsOut)
    END DO ! Row = 1, NumRowsOut, 1

    IF (ALLOCATED(Array2DOut)) DEALLOCATE(Array2DOut)

    CLOSE (UNIT = 10, STATUS = "KEEP", IOSTAT = CloseStat)
    IF (CloseStat /= 0) STOP "*** Cannot close file ***"

  END SUBROUTINE Write3DArrayAs2DArray



END MODULE ModArray
